% afni_dig_brik
% Christopher J. Whalen
% University of Illinois
% 07/01/04 - Date Created
% 08/09/04 - Last Modified
% 05/12/11 - fixed "tai_fid" to refer to "tai.xyz" and removed "tai_origin"
%            (no idea how this routine could ever have worked) - elm
%
% Purpose - to create a brik of the digitized points and the tai fiducials

function [dig_brik] = afni_dig_brik(choice,xyz,dig_brik,src,det,tai,opt_label,mri_dim)

fprintf('\nCreating Optode BRIK ...\n');
% define acpc.BRIK size and orientation

xyz = round(xyz);
xyz = xyz - 1;           % Subtracting 1 mm to go back into AFNI coords (?)
tai_fid = tai.xyz - 1;   % Subtracting 1 mm to go back into AFNI coords (?)

% SRC - Red
% color 1 voxel on all sides of dig pts of src (cubes of 27) 
for i_src=1:length(src)
    for i_opt_label=1:length(opt_label)
       if strcmp(src(i_src),opt_label(i_opt_label))==1
           x=xyz(i_opt_label,1);  
           y=xyz(i_opt_label,2);
           z=xyz(i_opt_label,3);
           if x<2;x=2;end;if y<2;y=2;end;if z<2;z=2;end
           if x>(mri_dim.x-2);x=mri_dim.x-2;end;if y>(mri_dim.y-2);y=mri_dim.y-2;end;if z>(mri_dim.z-2);z=mri_dim.z-2;end;
           dig_brik(x-1:x+1,y-1:y+1,z-1:z+1)=128;
       end
    end
end

% DET - Yellow
% color 1 voxel on all sides of dig pts of det (cubes of 27) 
for i_det=1:length(det)
    for i_opt_label=1:length(opt_label)
       if strcmp(det(i_det),opt_label(i_opt_label))==1
           x=xyz(i_opt_label,1);
           y=xyz(i_opt_label,2);
           z=xyz(i_opt_label,3);
           if x<2;x=2;end;if y<2;y=2;end;if z<2;z=2;end
           if x>(mri_dim.x-2);x=mri_dim.x-2;end;if y>(mri_dim.y-2);y=mri_dim.y-2;end;if z>(mri_dim.z-2);z=mri_dim.z-2;end;
           dig_brik(x-1:x+1,y-1:y+1,z-1:z+1)=88;
       end
    end
end

% XYZ fid - Red
% color 1 voxel on all sides of dig pts of xyz fid (cubes of 27) 
for i_fid=1:3
       x=xyz(i_fid,1);
       y=xyz(i_fid,2);
       z=xyz(i_fid,3);
       if x<2;x=2;end;if y<2;y=2;end;if z<2;z=2;end
       if x>(mri_dim.x-2);x=mri_dim.x-2;end;if y>(mri_dim.y-2);y=mri_dim.y-2;end;if z>(mri_dim.z-2);z=mri_dim.z-2;end;
       dig_brik(x-1:x+1,y-1:y+1,z-1:z+1)=128;
end


% XYZ Dummy - Light Blue (these will be removed before this subrouting starts if so desired)
% color 1 voxel on all sides of dig pts of xyz dummy (cubes of 27)
for i_opt_label=3:length(opt_label) % ignoring fiducials - they are not dummy points
    dum_pt_index = 1; % This will be 1 if the point is a dummy point and 0 if not dummy point
    for i_src=1:length(src)
        if strcmp(src(i_src),opt_label(i_opt_label))==1 % Seeing if point is src point and not a dummy point
            dum_pt_index = 0;
        end
    end
    for i_det=1:length(det)
        if strcmp(det(i_det),opt_label(i_opt_label))==1 % Seeing if point is det point and not a dummy point
            dum_pt_index = 0;
        end
    end
    if dum_pt_index == 1 % if dummy point color a diff color
       x=xyz(i_opt_label,1);
       y=xyz(i_opt_label,2);
       z=xyz(i_opt_label,3);
       if x<2;x=2;end;if y<2;y=2;end;if z<2;z=2;end
       if x>(mri_dim.x-2);x=mri_dim.x-2;end;if y>(mri_dim.y-2);y=mri_dim.y-2;end;if z>(mri_dim.z-2);z=mri_dim.z-2;end;
       dig_brik(x-1:x+1,y-1:y+1,z-1:z+1)=45;
    end
end

% TAI Fid - Green
% color 1 voxel on all sides of dig pts of tai fid (cubes of 27) (3 voxels
% for RPA (cubes of 343) 
for i_fid=1:3
       x=tai_fid(i_fid,1);
       y=tai_fid(i_fid,2);
       z=tai_fid(i_fid,3);
       if x<2;x=2;end;if y<2;y=2;end;if z<2;z=2;end
       if x>(mri_dim.x-2);x=mri_dim.x-2;end;if y>(mri_dim.y-2);y=mri_dim.y-2;end;if z>(mri_dim.z-2);z=mri_dim.z-2;end;
       if i_fid == 3
           if x<4;x=4;end;if y<4;y=4;end;if z<4;z=4;end
           if x>(mri_dim.x-4);x=mri_dim.x-4;end;if y>(mri_dim.y-4);y=mri_dim.y-4;end;if z>(mri_dim.z-4);z=mri_dim.z-4;end;
           dig_brik(x-3:x+3,y-3:y+3,z-3:z+3)=4;
       else    
           if x<2;x=2;end;if y<2;y=2;end;if z<2;z=2;end
           if x>(mri_dim.x-2);x=mri_dim.x-2;end;if y>(mri_dim.y-2);y=mri_dim.y-2;end;if z>(mri_dim.z-2);z=mri_dim.z-2;end;
           dig_brik(x-1:x+1,y-1:y+1,z-1:z+1)=4;
       end
end

% Remove Back of Head Plane
%
if choice.rm_pts == 2 | choice.rm_pts == 3
    range = max(xyz(:,2))-xyz(1,2);
    init_test_size = size(xyz,1);
    fprintf('\nEnter Back of Head Threshold Fraction:\n'); %(Threshold Fraction is the fraction of the distance from the most posterior elp point to the nasion elp point from which to remove posterior points.\n(threshold fraction of zero is keeping all points. threshold fraction of 1 is removal of all points. threshold of 0.30 (default) is removing the points lying in the posterior 30% of the scalp)\n'); 
    thresh_fraction = input('','s');
    if strcmp(thresh_fraction,'') ==1
      thresh_fraction = 0.30; % Removing the most posterior 3/10 pts
    else
      thresh_fraction = str2num(thresh_fraction);
    end
    thresh_plane_AP = xyz(1,2) + range*(1-thresh_fraction);
    y = round(thresh_plane_AP);
    for x=1:mri_dim.x
        for z=1:mri_dim.z
            dig_brik(x,y,z)=80;
        end
    end
end

% Remove Face Plane
%
if choice.rm_pts == 1 | choice.rm_pts == 3
    range = max(xyz(:,2))-xyz(1,2);
    init_test_size = size(xyz,1);
    fprintf('\nEnter Face Threshold Fraction:\n'); %(Threshold Fraction is the fraction of the distance from the most posterior elp point to the nasion elp point from which to remove posterior points.\n(threshold fraction of zero is keeping all points. threshold fraction of 1 is removal of all points. threshold of 0.30 (default) is removing the points lying in the posterior 30% of the scalp)\n'); 
    thresh_fraction = input('','s');
    if strcmp(thresh_fraction,'') ==1
      thresh_fraction = 0.10; % Removing the most posterior 3/10 pts
    else
      thresh_fraction = str2num(thresh_fraction);
    end
    thresh_plane_AP = xyz(1,2) + range*thresh_fraction;
    y = round(thresh_plane_AP);
    for x=1:mri_dim.x
        for z=1:mri_dim.z
            dig_brik(x,y,z)=80;
        end
    end
end

% % color 1 voxel on all sides of dig pts of det (cubes of 27) 
% for i_det=1:length(det)
%     for i_opt_label=1:length(opt_label)
%        if strcmp(det(i_det),opt_label(i_opt_label))==1
%            x=xyz(i_opt_label,1);
%            y=xyz(i_opt_label,2);
%            z=xyz(i_opt_label,3);
%            if x<2;x=2;end;if y<2;y=2;end;if z<2;z=2;end
%            if x>(mri_dim.x-2);x=mri_dim.x-2;end;if y>(mri_dim.y-2);y=mri_dim.y-2;end;if z>(mri_dim.z-2);z=mri_dim.z-2;end;
%            dig_brik(x-1:x+1,y-1:y+1,z-1:z+1)=64;
%        end
%     end
% end
